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Abstract 

We propose an optimal control framework for persistent monitoring problems where the objective 
is to control the movement of mobile nodes to minimize an uncertainty metric in a given mission 
space. For multi agent in a one-dimensional mission space, we show that the optimal solution is 
obtained in terms of a sequence of switching locations and waiting time on these switching points, 
thus reducing it to a parametric optimization problem. Using Infinitesimal Perturbation Analysis 
(IPA) we obtain a complete solution through a gradient-based algorithm. We also discuss a receding 
horizon controller which is capable of obtaining a near-optimal solution on-the-fly. 



1 Introduction 

Enabled by recent technological advances, the deployment of autonomous agents that can cooperatively 
perform complex tasks is rapidly becoming a reality. In particular, there has been considerable progress 
reported in the literature on robotics and sensor networks regarding coverage control |l]-[3j, surveil- 
lance plB) and environmental sampling (6j|7j missions. In this paper, we are interested in generating 
optimal control strategies for persistent monitoring tasks; these arise when agents must monitor a dy- 
namically changing environment which cannot be fully covered by a stationary team of available agents. 
Persistent monitoring differs from traditional coverage tasks due to the perpetual need to cover a chang- 
ing environment, i.e., all areas of the mission space must be visited infinitely often. The main challenge 
in designing control strategies in this case is in balancing the presence of agents in the changing en- 
vironment so that it is covered over time optimally (in some well-defined sense) while still satisfying 

*The authors' work is supported in part by NSF under Grant EFRI-0735974, by AFOSR under grant FA9550-09- 1-0095, 
by DOE under grant DE-FG52-06NA27490, by ONR under grant N000 14-09- 1-1051 and by ARO under grant W911NF-1 1-1- 
0227. 
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sensing and motion constraints. Examples of persistent monitoring missions include surveillance and 
theft prevention in a building, patrol missions with unmanned vehicles, and environmental applications 
where routine sampling of an area is involved. 

In this paper, we address the persistent monitoring problem by proposing an optimal control framework 
to drive agents so as to minimize a metric of uncertainty over the environment. In coverage control |2j[3j, 
it is common to model knowledge of the environment as a non-negative density function defined over the 
mission space, and usually assumed to be fixed over time. However, since persistent monitoring tasks 
involve dynamically changing environments, it is natural to extend it to a function of both space and time 
to model uncertainty in the environment. We assume that uncertainty at a point grows in time if it is not 
covered by any agent sensors. To model sensor coverage, we define a probability of detecting events at 
each point of the mission space by agent sensors. Thus, the uncertainty of the environment decreases 
with a rate proportional to the event detection probability, i.e., the higher the sensing effectiveness is, 
the faster the uncertainty is reduced.. 

While it is desirable to track the value of uncertainty over all points in the environment, this is generally 
infeasible due to computational complexity and memory constraints. Motivated by polling models in 
queueing theory, e.g., spatial queueing j8), j9J, and by stochastic flow models 1 10 1, we assign sampling 



points of the environment to be monitored persistently (this is equivalent to partitioning the environment 
into a discrete set of regions.) We associate to these points "uncertainty queues" which are visited by 
one or more "servers". The growth in uncertainty at a sampling point can then be viewed as a flow into 
a queue, and the reduction in uncertainty (when covered by an agent) can be viewed as the queue being 
visited by mobile servers as in a polling system. Moreover, the service flow rates depend on the distance 
of the sampling point to nearby agents. From this point of view, we aim to control the movement of the 
servers (agents) so that the total accumulated "uncertainty queue" content is minimized. 

Control and motion planning for agents performing persistent monitoring tasks have been studied in the 
literature. In pi the focus is on sweep coverage problems, where agents are controlled to sweep an area. 
In j6j[TTJ a similar metric of uncertainty is used to model knowledge of a dynamic environment. In |TTJ, 
the sampling points in a 1 -dimensional environment are denoted as cells, and the optimal control policy 
for a two-cell problem is given. Problems with more than two cells are addressed by a heuristic policy. 
In [6], the authors proposed a stabilizing speed controller for a single agent so that the accumulated 
uncertainty over a given path in the environment is bounded, along with an optimal controller that min- 
imizes the maximum steady-state uncertainty, assuming that the agent travels along a closed path and 
does not change direction. The persistent monitoring problem is also related to robot patrol problems, 
where a team of robots are required to visit points in the workspace with frequency constraints 1 12-[T4|. 



Our ultimate goal is to optimally control a team of cooperating agents in a 2 or 3-dimensional environ- 
ment. The contribution of this paper is to take a first step toward this goal by formulating and solving 
an optimal control problem for a team of agents moving in a 1 -dimensional mission space described by 
an interval [0,L] C K in which we minimize the accumulated uncertainty over a given time horizon and 
over an arbitrary number of sampling points. Even in this simple case, determining a complete explicit 
solution is computationally hard. However, we show that the problem can be reduced to a parametric 
optimization problem. In particular, the optimal trajectory of each agent is to move at full speed until 
it reaches some switching point, dwell on the switching point for some time (possibly zero), and then 
switch directions. In addition, we prove that all agents should never reach the end points of the mis- 
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sion space [0,L]. Thus, each agent's optimal trajectory is fully described by a set of switching points 
{8\, . . . , Qk} and associated waiting times at these points, {w\,. . . ,wk}- As a result, we show that the 
behavior of the agents operating under optimal control is described by a hybrid system. This allows 



us to make use of generalized Infinitesimal Perturbation Analysis (IPA), as presented in 1 15 1, 1 16 1, to 
determine gradients of the objective function with respect to these parameters and subsequently obtain 
optimal switching locations and waiting times that fully characterize an optimal solution. It also allows 
us to exploit robustness properties of IPA to extend this solution approach to a stochastic uncertainty 
model. Our analysis establishes the basis for extending this approach to a 2-dimensional mission space 
(in ongoing research). In a broader context, our approach brings together optimal control, hybrid sys- 
tems, and perturbation analysis techniques in solving a class of problems which, under optimal control, 
can be shown to behave like hybrid systems characterized by a set of parameters whose optimal values 
deliver a complete optimal control solution. 

The rest of the paper is organized as follows. Section [2] formulates the optimal control problem. Section 
[3] characterizes the solution of the optimal control problem in terms of two parameter vectors specifying 
switching points in the mission space and associated dwelling times at them. Using IPA in conjunction 
with a gradient-based algorithm, a complete solution is also provided. Section|4]provides some numerical 
results and Section |5]concludes the paper. 



2 Persistent Monitoring Problem Formulation 

We consider N mobile agents moving in a 1 -dimensional mission space of length L, for simplicity taken 
to be an interval [0,L] Cl. Let the position of the agents at time t be s n (t) € [0,L],n= 1, ... ,N, following 
the dynamics: 

s n (t) = u n (t) (1) 

i.e., we assume that the agent can control its direction and speed. Without loss of generality, after 
some rescaling with the size of the mission space L, we further assume that the speed is constrained by 
\u n — 1. n = 1 ? ■ • ■ ,N. For the sake of generality, we include the additional constraint: 

a<s(t) <b, a>0, b <L (2) 

over all t to allow for mission spaces where the agents may not reach the end points of [0,L], possibly due 
to the presence of obstacles. We also point out that the agent dynamics in ([I]) can be replaced by a more 
general model of the form s n (t) = g n (s n ) + b n u n {t) without affecting the main results of our analysis (see 



also Remark lin Section 3.1 ) Finally, an additional constraint may be imposed if we assume that the 
agents are initially located so that s„ (0) < s n+ \ (0), n = I, . . . ,N — I, and we wish to prevent them from 
subsequently crossing each other over all t: 

Sn(t)-s n +i(t)<0 (3) 

We associate with every point x £ [0,L] a function p n (x,s„) that measures the probability that an event 
at location x is detected by agent n. We also assume that p n (x,s n ) = 1 if x = s n , and that p n (x,s n ) 
is monotonically nonincreasing in the distance \x — s n \ between x and s n , thus capturing the reduced 
effectiveness of a sensor over its range which we consider to be finite and denoted by r n (this is the same 
as the concept of "sensor footprint" found in the robotics literature.) Therefore, we set p n (x,s n ) = 
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when \x — s n \ > r n . Although our analysis is not affected by the precise sensing model p n (x, s n ), we will 
limit ourselves to a linear decay model as follows: 

Pn(x,S n ) = \ ^- Sn \ (4) 

[ 0, if\x-s n \ >r n 

Next, consider a set of points {a,}, i = 1, . . . ,M, a, G [0,L], and associate a time- varying measure of 
uncertainty with each point a,, which we denote by Ri(t). Without loss of generality, we assume < 
CL\ < ■■ ■ < GCm < L and, to simplify notation, we set p n ,i(s n (t)) = p n ((Xi,s n (t)). This set may be selected 
to contain points of interest in the environment, or sampled points from the mission space. Alternatively, 
we may consider a partition of [0,L] into M intervals whose center points are a,- = 2M > 1 ' = lj • • • 
We can then set p n (x,s n (?)) = p n j(s n (/)) for all x G [a,- — + ^]. Therefore, the joint probability 
of detecting an event at location x G [a,- — ^,ce,- + j^] by all the AT agents simultaneously (assuming 
detection independence) is: 

Q 

PiKt)) = l-H[l-Pn,i(Sn(t))} (5) 

n=l 

where we set s(t) = [si (t) , ... ,s^(t)] T . We define uncertainty functions Ri(t) associated with the in- 
tervals [a, ■ — ^ , a,- + 2^], i = 1, . . . ,M, so that they have the following properties: (i) Ri(t) increases 
with a prespecified rate A, if Pi (s(?)) = 0, (//) Ri(t) decreases with a fixed rate B if Z 5 ,- (s(f)) = 1 and 
(Hi) Rt(t) > for all t. It is then natural to model uncertainty so that its decrease is proportional to the 
probability of detection. In particular, we model the dynamics of /?,(?), i = 1, . . . ,M, as follows: 

HR i (t)=0,A i <BFl(s(t)) 
l{ ) \ Ai-BPi(s(t)) otherwise K ' 

where we assume that initial conditions /?/(0), i = 1, . . . ,M, are given and that B > A,- > (thus, the 
uncertainty strictly decreases when there is perfect sensing Pi (s(?)) = 1.) 

Viewing persistent monitoring as a polling system, each point a,- (equivalently, ith interval in [0,L]) is 
associated with a "virtual queue" where uncertainty accumulates with inflow rate A,-. The service rate 
of this queue is time- varying and given by BPi (s(t)), controllable through the agent position at time t. 
Figure [T] illustrates this polling system when N = \. This interpretation is convenient for characterizing 
the stability of such a system over a mission time T: For each queue, we may require that fj Ai < 
fo Bpi(s(t))dt. Alternatively, we may require that each queue becomes empty at least once over [0,T]. 
We may also impose conditions such as Ri(T) < 7? max for each queue as additional constraints for our 
problem so as to provide bounded uncertainty guarantees, although we will not do so in this paper. Note 
that this analogy readily extends to 2 or 3-dimensional settings. 

The goal of the optimal persistent monitoring problem we consider is to control the movement of the ,/V 
agents through u„ (t) in ([TJ so that the cumulative uncertainty over all sensing points {a,}, i = 1, . . . ,M 
is minimized over a fixed time horizon T. Thus, setting u (t) = [u\ (t) , . . . ,un (t)] we aim to solve the 
following optimal control problem PI: 

l r T M 

u(f) T Jo f-[ 

subject to the agent dynamics ([T]), uncertainty dynamics ([6]), control constraint < 1, t G [0,T], 

and state constraints (jl]), t G [0,T]. Note that we require a < r n and b > L — r m , for at least some 
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Figure 1 : A queueing system analog of the persistent monitoring problem. 



n,m= 1, . . . ,N; this is to ensure that there are no points in [0,L] which can never be sensed, i.e., any i 
such that a, < a — r n or a,- > b + r n would always lie outside any agent's sensing range. We will omit the 
additional constraint Q from our initial analysis, but we will show that, when it is included, the optimal 
solution never allows it to be active. 



3 Optimal Control Solution 



3.1 Hamiltonian analysis 



We first characterize the optimal control solution of problem PI and show that it can be reduced to a 
parametric optimization problem. This allows us to utilize an Infinitesimal Perturbation Analysis (IPA) 
gradient estimation approach fT"5| to find a complete optimal solution through a gradient-based algorithm. 
We define the state vector x (?) = [s\ (?) , . . . , tfjv(0 (0 > • ■ • (0] T an d the associated costate vector 
X (?) = [A Vl (?),... > A JA , (t) , Ai (t) , . . . , %m (0] T - m view of the discontinuity in the dynamics of Rt(t) in 
([6]>, the optimal state trajectory may contain a boundary arc when Ri(t) = for any i; otherwise, the 
state evolves in an interior arc. We first analyze the system operating in such an interior arc and omit the 
constraint Q as well. Using ([T]) and ([6]), the Hamiltonian is 

M N M 

H (x, X ,u) = + £ K (0 «n (t) + £ Xi (t)Ri(t) (8) 

i=l n=\ i=\ 

and the costate equations A = — are 

1, i = l,...,M (9) 

- I Xi(t)Y\[i-p d Mi{m+- I mWy-p*MW>)\ m 

rn ieF7At) d + n r,1 ieFi(t) d+n 



K (0 



dH 

dRi(t) 

dH 
ds„ (t) 
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where we have used ([4]), and the sets F n (?) and F„(t) are defined as 

Fn(t) = {i : s„ (t) - r„ < a t < s„ (?)} (11) 
F+O) = {i : s„ (?) < a,- < s„ (?) + r„} 

for « = 1, . . . ,iV. Note that F„{t), F,j~(?) identify all points 05; to the left and right of s n (?) respectively 
that are within agent ti's sensing range. Since we impose no terminal state constraints, the boundary 
conditions are A,- (T) = 0, i = 1, . . . ,M and X Sn (T) = 0, n = l,...,N. Applying the Pontryagin minimum 
principle to ^ with u*(?), ? G [0, T), denoting an optimal control, we have 

H(x*,X*,u*)= min H(x,k,u) 

K„e[-l,l], n=L...,N 

and it is immediately obvious that it is necessary for an optimal control to satisfy: 

u*(t) = { 1 ifA -W <0 (12) 
U » {t) \ -1 ifA i „(0>0 {U) 

This condition excludes the possibility that A iS -„ (?) = over some finite singular intervals (17) . We will 
show that if s n (?) = a > or s n (t) = b < L, then X Sji (t) = for some n € {1, . . . , A^} may in fact exist for 
some finite arc; otherwise X Sn (t) = can arise only when u n (t) = 0. 

The implication of ^ with A,- (J) = is that A,- (t) = T -t for all f G [0, 7 1 ] and all z = 1,...,M and that 
A, (?) is monotonically decreasing starting with A, (0) = T . However, this is only true if the entire optimal 
trajectory is an interior arc, i.e., all Rt{t) > constraints for all i = 1 , . . . ,M remain inactive. On the other 



hand, looking at ( 10 ), observe that when the two end points, and L, are not within the range of an agent, 
we have |F W ~(?)| = |F„ + (?)|, since the number of indices i satisfying s n (?) — r n < a,- < s n (?) is the same 
as that satisfying s n (?) < a, < s n (?) + r„. Consequently, for the one-agent case N =1, ( flo] ) becomes 

k(t) = -- i m+- i m d3) 

1 '-effW 1 teF+« 



and X S{ (?) = since the two terms in ( 13 1 will cancel out, i.e., A. Vl (?) remains constant as long as this 
condition is satisfied and, in addition, none of the state constraints /?;(?) > 0, i = 1, . .. ,M, is active. Thus, 
for the one agent case, as long as the optimal trajectory is an interior arc and X Sl (?) < 0, the agent moves 
at maximal speed u\ (?) = 1 in the positive direction towards the point si = b. If X Sl (?) switches sign 
before any of the state constraints /?,(?) > 0, i = 1, . . . ,M, becomes active or the agent reaches the end 
point s\ = b, then u\ (?) = — 1 and the agent reverses its direction or, possibly, comes to rest. 

In what follows, we examine the effect of the state constraints which significantly complicates the analy- 
sis, leading to a challenging two-point-boundary- value problem. However, we will establish the fact that 
the complete solution boils down to determining a set of switching locations over [a, b] and waiting times 
at these switching points, with the end points, and L, being always infeasible on an optimal trajectory. 
This is a much simpler problem that we are subsequently able to solve. 

We begin by recalling that the dynamics in ^ indicate a discontinuity arising when the condition /?, (?) = 
is satisfied while /?,(/) = A, — BPi (s(?)) < for some i = 1, . . . ,M. Thus, = defines an interior 
boundary condition which is not an explicit function of time. Following standard optimal control analysis 
[ 17 ], if this condition is satisfied at time ? for some j G {1, . . . ,M}, 

//(x(?-),A(?-),u(?-)) =//(x(?+),A(? + ),u(?+)) (14) 
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where we note that one can choose to set the Hamiltonian to be continuous at the entry point of a boundary 
arc or at the exit point. Using ([8]) and ([6]), ( [14] ) implies: 

f a; (O K (r) + a; (?-) [Aj (?) - BPMm = £ a* (?+) «; (?+) (is) 

In addition, A* (f~) = X* (t+) for all n = 1,. . . ,/V and A* (f~) = A* (?+) for all j / j, but Xj (?) may 
experience a discontinuity so that: 

a;(?-)=a;(? + )-% (i6) 

where TCj > is a multiplier associated with the constraint —Rj(t) < 0. Recalling (12 1, since A* (?) 
remains unaffected, so does the optimal control, i.e., u*(t~) = u*(t + ). Moreover, since this is an entry 



point of a boundary arc, it follows from (|6j) that Aj — BPj (s(t)) < 0. Therefore, ( 15 I and ( 16 1 imply that 

A;(r)=0, X}{t + )=xj>0. 

Thus, A,- (?) always decreases with constant rate — 1 until /?, (?) = is active, at which point A, (?) jumps 
to a non-negative value 71, and decreases with rate — 1 again. The value of 7T; is determined by how long 
it takes for the agents to reduce /?, (?) to once again. Obviously, 

Aj (?) > 0, / = 1, . . . ,M, ? G [0, T] (17) 

with equality holding only if ? = T, or? = t with R, (to) = 0, /?;(/') > 0, where?' G [?o-5,? ), 5 > 0. The 



actual evaluation of the costate vector over the interval [0, T] requires solving ( 10 ), which in turn involves 
the determination of all points where the state variables R,-(t) reach their minimum feasible values /?,•(?) = 
0, i = 1, . . . ,M. This generally involves the solution of a two-point-boundary- value problem. However, 
our analysis thus far has already established the structure of the optimal control ( fp2| ) which we have seen 
to remain unaffected by the presence of boundary arcs when /?;(?) = for one or more i = 1, . . . ,M. We 
will next prove some additional structural properties of an optimal trajectory, based on which we show 
that it is fully characterized by a set of non-negative scalar parameters. Determining the values of these 
parameters is a much simpler problem that does not require the solution of a two-point-boundary-value 
problem. 

Let us turn our attention to the constraints s„ (?) > a and s„ (?) < b and consider first the case where a = 0, 
b = L, i.e., the agents can move over the entire [0,L]. We shall make use of the following technical 
condition: 

Assumption 1: For any n = I,... ,N, i= 1, . . . ,M, t G (0, T), and any e > 0, if s„(t) = 0, s n (t — e) > 0, 
then either Ri(x) > for all x G [?-£,?] or Ri(x) = for all X G [?-£,?]; if s„(t) =L,s n (t — e) < L,then 
either Rj(x) > for all X G [? — £,?] or /?,(t) = for all X G [?- £,?]. 

This condition excludes the case where an agent reaches an endpoint of the mission space at the exact 
same time that any one of the uncertainty functions reaches its minimal value of zero. Then, the following 
proposition asserts that neither of the constraints s n (t) > and s n (t) <L can become active on an optimal 
trajectory. 

Proposition 3.1 Under Assumption 1, if a = 0, b = L, then on an optimal trajectory: s*(t) ^ and 
s*(t)^Lforallte (0,T), «€{!,... ,N}. 
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Proof. Suppose at t = to < T an agent reaches the left endpoint, i.e., s* (to) = 0, s* n Uq) > 0. We will then 
establish a contradiction. Thus, assuming s* (to) = 0, we first show that X* Uq) = by a contradiction 
argument. Assume that X* (?q~) / 0, in which case, since the agent is moving toward s„ = 0, we have 



u* \ t Q ) = — 1 and X* \t Q ) > from ( 12). Then, X* (t) may experience a discontinuity so that 

K( t o)=KA t O + )-*n (18) 

where % n ^ is a scalar constant. It follows that A* (?^) = A* (t^) + 7t n > 0. Since the constraint 
s n (?) = is not an explicit function of time, we have 

K^)<{h)=KM)<{4) d9) 

On the other hand, u* (Jq) ^ 0, since agent n must either come to rest or reverse its motion at s n = 0, 



hence A* Uq) u* n (t^) ^ 0. This violates ( 19 1, since A* (t n ) u* (t Q ) < 0. This contradiction implies that 



A s * (t Q ) = 0. Next, consider (flOj) and observe that in ( 1 1 ) we have F n (to) = 0, since a,- > s* (to) = for 



all i = 1, . . . ,M. Therefore, recalling ([17]), it follows from ([T0J) that 

M*o) = r I ^(Ollt 1 -^^))] ^° 

r "/er,l(r -) ^ 

Under Assumption 1, there exists 5i > such that during the interval (to — 8\,to) no Rj (t) > becomes 
active, hence no Aj(f) encounters a jump for / = 1, . . . ,M. It follows that A ( *(f) > for i G F n + (f) and 
A* (f) is continuous with A 4 * (?) > for t G (?o — 5i,fo). Again, since s* (to) = 0, there exists some 
&2 < 8\ such that for t G (to — 82, t ), we have u* n (t) < and A* (t) > 0. Thus, for t G (to — &i, to), we 
have A* (f) > and A* ; (f) > 0. This contradicts the fact we already established that A s * (f^) =0 and 
we conclude that s* (t) 7^ for all t G [0,T], n = 1, . . . ,JV. Using a similar line of argument, we can also 
show that s* (t) / L. ■ 

Proposition 3.2 Tfa > and (or) b < L, then on an optimal trajectory there exist finite length intervals 
[to,t\] such that s n (t) = a and (or) s n (t) = b,for some n G {!,. . . ,A^}, t G [fo; f i]> ^ k) <h <T. 



Proof. Proceeding as in the proof of Proposition 3.1 when s* n (to) = a we can establish ( 19 ) and the fact 



that A* (? ) = 0. On the other hand, u* (t^) ^ 0, since the agent must either come to rest or reverse its 
motion at s n (to) = a. In other words, when s n (to) = a on an optimal trajectory, (fl9]) is satisfied either 



with the agent reversing its direction immediately (in which case t\ = to and A* (t^) = 0) or staying on 
the boundary arc for a finite time interval (in which case t\ > to and u* n (t) = for t G [^o 5 f 1 ] ) - T ne exact 
same argument can be applied to s n (t) = b. ■ 

The next result establishes the fact that on an optimal trajectory, every agent either moves at full speed 
or is at rest. 

Proposition 3.3 On an optimal trajectory, either u* (t) = ±1 ifX* n (t) 7^ 0, or u* n (t) = if X* (t) = Ofor 
te[0,T],n = l,...,N. 



Proof. When A* (t) ^ 0, we have shown in ( 12 1 that u* n (t) = ±1, depending on the sign of A* (t). Thus, 



it remains to consider the case A* (t) = for some t G [?i,?2], where < t\ < ?2 < T. Since the state is 
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in a singular arc, A* (?) does not provide information about u* (?). On the other hand, the Hamiltonian in 
^ is not a explicit function of time, therefore, setting H (x*,A*,u*) = H* , we have^- = 0, which gives 



N 



M 



M 



AV* M N 

E**(o+E w «: (o+E ^: (o «: w + e v o # w + i a * w *? w = o (20) 



d? 



1=1 



n=l 



i=l 



Define 5(f) = {«|A Jn (?) = 0,n = 1, . . . ,Af} as the set of indices of agents that are in a singular arc and 
S(?) = {n\X Sn (?) ^ 0,n = 1, . . . ,iV} as the set of indices of all other agents. Thus, A* (?) = 0, A* (?) = 
for ? G [?i ,?2] ,n € 5 (?). In addition, agents move with constant full speed, either 1 or — 1 , so that u* (?) = 0, 
n G S(t). Then, (20 1 becomes 



E[i+V «]*?(*)+ E (0«;(0+Ev (0*?(0=o 



d? 



(21) 



1=1 



1=1 



From (|9J, A* (?) = — 1, i = 1,... ,M, so 1 + A* (f) = 0, leaving only the last two terms above. Note that 

K, ( t ) = -M) and writing %® = dJ tr we § et: 



E<«^ + E VW^^o 

neS(t) i=l,i?,^0 



Recall from (6) that when /?,•(?) / Owe have /?,■(?) = A,-B[l - ]^[ [1 -p,-(j„ (?))]], so that 



n=l 



—---B 2. a H(i -MM'))) 

^nlU i=l,R,^0 ff! »l'J d^n 



dt 



„=1 W rf^„ 



which results in 

M 

s E vw 

i= 1,R,^0 



d Pi (s* n (t)) 



£ " : « n (! - « w (0)) - e < it) n a - * w (0)) 

„=1 °^iAU ^„ 



nSS(f) 



d Pi (s*(t)) 



-b £ VW E <W^Pii(i-ftWW)) = o 

i=l,R,.^0 neS(t) as n\ l ) d+n 



(22) 



Note that -^^m^ = ±^ or 0, depending on the relative position of s\ (?) with respect to a,-. Moreover, 
(22 1 is invariant to M or the precise way in which the mission space [0,L] is partitioned, which implies 
that 

V(0 E <(O^M^n(i-M4(0)) = o 

for all i = 1, ... ,M, f G [?i,?2] • Since A* (?) = — 1, i = 1,... ,M, it is clear that to satisfy this equality we 
must have u* (?) = for all ? G [?i,?2] €S(t). In conclusion, in a singular arc with A* (0=0 f° r some 
n G { 1 , . . . , N} , the optimal control is u* (?) = 0. ■ 



9 



Next, we consider the case where the additional state constraint ([3]) is included. We can then prove that 
this constraint is never active on an optimal trajectory, i.e., agents reverse their direction before making 
contact with any other agent. 

Proposition 3.4 If the constraint Q is included in problem PI, then on an optimal trajectory, s*(t) 7^ 
s* +1 (t)farte(0,T},n = l,...,N-l. 



Proof. Suppose at t = to < T we have s* (to) = s* +l (to), for some n = 1, . . . , N— 1. We will then establish 
a contradiction. First assuming that both agents are moving (as opposed to one being at rest) toward each 



other, we have u* n {t^ \ = 1 and u* +l (t^ = —1. From (12i and Prop 3.3 we know A s * (f^) < and 
A* (?q~) > 0. When the constraint s n (?) — s n+ i (t) < is active, A* (t) and A* (f^) may experience a 
discontinuity so that 

K+ 1 fo)=*L. l (4)-* (24) 

where % ^ is a scalar constant. It follows that A* (t^) = A* (t Q ) —K<0 and X* (t$) = A* +i (t Q ) + 
% > 0. Since the constraint s n (t) — s n+ \ (t) < is not an explicit function of time, we have 

K ( { o) < {to) +K +1 ih) < + i ih) = K {*o) < {*o) +K, +i {to) ('o + ) (25) 

On the other hand, u* (fj") ^ and u* n+l (t^) ^ 0, since agents n and n + 1 must either come to rest or 
reverse their motion after making contact, hence X* n (f^) u* (f^) + A^* (^) u* n+l {t£) ^ 0. This violates 



(25 1, since A 4 * (t Q ) u* (t Q ) + X* (t^) u* n+l (t^) < 0. This contradiction implies that s n (t) —s n+ \ (t) = 
cannot be active and we conclude that s* n (t) ^ s* n+l (t) for t £ [0,T], n = 1, . . . ,Af — 1. Moreover, if one 
of the two agents is at rest when s* (to) = s* n+l (to), the same argument still holds since it is still true that 

K ^0) < ^0) +K +l {t Q ) K+i ^0) < o- ■ 

Based on this analysis, the optimal control u* (t) depends entirely on the sign of X* (t) and, in light 



of Propositions 3.1 3.3 the solution of the problem reduces to determining: (/) switching points in 
[0,L] where an agent switches from u* (t) = ±1 to either =pl or 0; or from u* (t) = to either ±1, and 
(ii) if an agent switches from u* (t) = ±1 to 0, waiting times until the agent switches back to a speed 
u* (t) = ±1. In other words, the full solution is characterized by two parameter vectors for each agent 
n: d n = . . , Qn,r„] T and w„ = \w n ,\ ■ ■ ■ ,Wn,r n ] T , where 6 n * G (0,L) denotes the <^th location where 

agent n changes its speed from ±1 to and w n p > denotes the time (which is possibly null) that 
agent n dwells on 8 n Note that Y n is generally not known a priori and depends on the time horizon 
T. In addition, we always assume that agent n reverses its velocity direction after leaving the switching 
point 6 n p with respect to the one it had when reaching d n p . This seemingly excludes the possibility 
of an agent's control following a sequence 1,0,1 or —1,0,-1. However, these two motion behaviors 
can be captured as two adjacent switching points approaching each other: when \ n p — n t+\ \ ~^ 0> tne 
agent control follows the sequence 1,0, 1 or —1,0, —1, and the waiting time associated with u* (?) = is 



For simplicity, we will assume that s n (0) = 0, so that it follows from Proposition 3.1 that w*(0) = 1, 



n = l,..,,N. Therefore, n \ corresponds to the optimal control switching from 1 to 0. Furthermore, 
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Q n £ with £ odd (even) always corresponds to u*(t) switching from 1 to (—1 to 0.) Thus, we have the 
following constraints on the switching locations for all E, = 2, . . . ,F„: 



< n,S-i, if £ is even 
> d n,^-i, if k isodd - 

It is now clear that the behavior of each agent under the optimal control policy is that of a hybrid system 
whose dynamics undergo switches when w* (?) changes from ±1 to and from to =pl or when R{(t) 
reaches or leaves the boundary value /?,• = 0. As a result, we are faced with a parametric optimization 
problem for a system with hybrid dynamics. This is a setting where one can apply the generalized 



theory of Infinitesimal Perturbation Analysis (IPA) in [ 15 1, 1 16] to conveniently obtain the gradient of 
the objective function / in ([7]) with respect to the vectors 6 and w, and therefore, determine (generally, 
locally) optimal vectors 6* and w* through a gradient-based optimization approach. Note that this is 
done on line, i.e., the gradient is evaluated by observing a trajectory with given 6 and w over [0, T] based 
on which 6 and w are adjusted until convergence is attained using standard gradient-based algorithms. 

Remark 1. If the agent dynamics in ([I]) are replaced by a model such as s n (t) = g n {sn) + b n u n (t), observe 



that ( 12b still holds. The difference lies in ( 10 1 which would involve a dependence on d8 ^f^ and further 
complicate the associated two-point-boundary-value problem. However, since the optimal solution is 
also defined by a parameter vectors d n = , . ■ • , 6n,r„] T an d w n = [w n ,i Wn,rJ T f° r eacn agent n, we 
can still apply the IPA approach presented in the next section. 



3.2 Infinitesimal Perturbation Analysis (IPA) 

Our analysis thus far has shown that, on an optimal trajectory, the agent moves at full speed, dwells on a 
switching point (possibly for zero time) and never reaches either boundary point, i.e., < s*(t) < L. Thus, 
the Tith agent's movement can be parameterized through Q n = [Q n \ 0«,r„] T an d w n = • • • , Wn,r„] T 
where Q n * is the <^th control switching point and w n p is the waiting time for this agent at the <i;th 
switching point. Therefore, the solution of problem PI reduces to the determination of optimal parameter 
vectors 6* and w*, n = 1 , . . . , N. As we pointed out, the agent's optimal behavior defines a hybrid system, 
and the switching locations translate to switching times between particular modes of this system. This is 



similar to switching-time optimization problems, e.g., [ |T8[ , |T9J, [20], except that we can only control 
a subset of mode switching times. We make use of IPA in part to exploit robustness properties that the 
resulting gradients possess pTJ; specifically, we will show that they do not depend on the uncertainty 
model parameters A,, i = 1,...,M, and may therefore be used without any detailed knowledge of how 
uncertainty affects the mission space. 



3.2.1 One agent solution with a = and b = L 

To maintain some notational simplicity, we begin with a single agent who can move on the entire mission 
space [0,L] and will then provide the natural extension to multiple agents and a mission space limited to 
[a,b] C [0,L]. We present the associated hybrid automaton model for this single-agent system operating 
on an optimal trajectory. Our goal is to determine VJ(d,w), the gradient of the objective function J in (jT]) 
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[/? 2 <i<a l +r,i^ >0 



a, <s<f) v R l >0 



[$<*«a;.,4 >0] 



i^(f)=X^(f)=i 

[a 1 -r<s<^,^ >0] 




i^(f) = / 2; i(f) = -l 
[ft t KS<a l +r,R 1 >0] 



V 2 

*=A „ 




11 

,S = <7 
! ► 
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1 


[a l <s<P 1 ,R 1 >0] 




[/? 1 <s<a l ,i^ >0] 





4©=j;,i(0=-i 

u-r<s </)..!> >0] 



where 



ft 



Figure 2: Hybrid automaton for each a,. Red arrows represent events when the control switches between 
1 and —1. Blue arrows represent events when Rj becomes 0. Black arrows represent all other events. 



with respect to 6 and w, which can then be used in a gradient-based algorithm to obtain optimal parameter 
vectors 6* and w*, n = l,...,N. We will apply IPA, which provides a formal way to obtain state and 
event time derivatives with respect to parameters of hybrid systems, from which we can subsequently 
obtaining VJ(d,w). 



Hybrid automaton model. We use a standard definition of a hybrid automaton (e.g., see |22|) as the 
formalism to model the system described above. Thus, let q G Q (a countable set) denote the discrete 
state (or mode) and x G X C W 1 denote the continuous state. Let v G T (a countable set) denote a 
discrete control input and u G U C R m a continuous control input. Similarly, let 8 G A (a countable set) 
denote a discrete disturbance input and d G D C W a continuous disturbance input. The state evolution 
is determined by means of (/) a vector field / : Q x X x U xD->X, (ii) an invariant (or domain) set 
Inv : Q x T x A — > 2 X , (Hi) a guard set Guard : gxgxTxA- > 2 X , and (iv) a reset function r : 
QxQxXxTxA— >X. The system remains at a discrete state q as long as the continuous (time-driven) 
state x does not leave the set Inv(q, v,8). Ifx reaches a set Guard (q, q',v,8) for some q' G Q, a discrete 
transition can take place. If this transition does take place, the state instantaneously resets to (q' ,x') 
where x' is determined by the reset map r(q,q' \x,v,8). Changes in v and 8 are discrete events that 
either enable a transition from q to q' by making sure x G Guard(q,q' ,v,8) or force a transition out of 
q by making sure x ^ Inv(q,v,8). We will classify all events that cause discrete state transitions in a 
manner that suits the purposes of IPA. Since our problem is set in a deterministic framework, 8 and d 
will not be used. 
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We show in Fig. [2] a partial hybrid automaton model of the single-agent system where a = and b = L. 
Since there is only one agent, we set s(t) = s\ (t), u(t) = u\ (?) and 6 = d\ for simplicity. Due to 
the size of the overall model, Fig. [2] is limited to the behavior of the agent with respect to a single 
a,-, i G {1, . . . ,M} and ignores modes where the agent dwells on the switching points (these, however, 



are included in our extended analysis in Section 3.2.2 ) The model consists of 14 discrete states (modes) 



and is symmetric in the sense that states 1—7 correspond to the agent operating with u(t) = 1, and states 
8 — 14 correspond to the agent operating with u(t) = — 1. States where u{t) = are omitted since we 
do not include the waiting time parameter w = w\ here. The events that cause state transitions can be 
placed in three categories: (i) The value of R,-(t) becomes and triggers a switch in the dynamics of ([6]). 
This can only happen when Ri(t) > and /?/(?) = A,- — Bpi(s(t)) < (e.g., in states 3 and 4), causing 
a transition to state 7 in which the invariant condition is Ri(t) = 0. (ii) The agent reaches a switching 
location, indicated by the guard condition s(t) = 6^ for any £ = 1, . . . ,T. In these cases, a transition 

results from a state z to z + 7 if z = 1, ,6 and to z — 7 otherwise. (Hi) The agent position reaches one of 

several critical values that affect the dynamics of Ri(t) while Ri(t) > 0. Specifically, when s(t) = a,- — r, 
the value of pi(s(t)) becomes strictly positive and Ri(t) = A; — Bpi(s(t)) > 0, as in the transition 1 — >2. 
Subsequently, when s(t) = a, — r(l — Aj/B), as in the transition 2—^3, the value of pi(s(t)) becomes 
sufficiently large to cause R,-(t) = Aj — Bpj(s(t)) < so that a transition due to /?,(?) = becomes feasible 
at this state. Similar transitions occur when s(t) = GC,, s(t) = a,- + r(l —Aj/B), and s(t) = a\ + r. The 
latter results in state 6 where Ri(t) = A,- > and the only feasible event is s(t) = 6^ , % odd, when a switch 
must occur and a transition to state 13 takes place (similarly for state 8). 

IPA review. Before proceeding, we provide a brief review of the IPA framework for general stochastic 
hybrid systems as presented in fI3] ]. The purpose of IPA is to study the behavior of a hybrid system state 
as a function of a parameter vector £ for a given compact, convex set 0Ct'. Let {Zk(6)}, k = 
l,...,K, denote the occurrence times of all events in the state trajectory. For convenience, we set To = 
and % + i = T. Over an interval [Tk(d),Zi c +i(d)), the system is at some mode during which the time- 
driven state satisfies x = fk(x, 6,t). An event at z^ is classified as (i) Exogenous if it causes a discrete 
state transition independent of 6 and satisfies = 0; (ii) Endogenous, if there exists a continuously 
differentiable function : MP x — > E such that z^ = min{? > T^-i : g^ (x(d,t) , d) = 0}; and (Hi) 
Induced if it is triggered by the occurrence of another event at time z m < z^. IPA specifies how changes 
in 6 influence the state x(6,t) and the event times Zk(Q) and, ultimately, how they influence interesting 
performance metrics which are generally expressed in terms of these variables. 

Given B = [6\,..., 0r] T , we use the Jacobian matrix notation: x'(t) = ^gl^ , z' k = , k = 1, . . . ,K, 
for all state and event time derivatives. It is shown in |15[ that x'(t) satisfies: 

d i , s dfk(t) i , ^ dfk(t) 

for t G [Zk,Zk+i) with boundary condition: 

x'(z+) = x'(z k -)+[f k _ 1 (z k -)-f k (z+)) z' k (28) 



for k = 0, . . . ,K. In addition, in (28 1, the gradient vector for each Zk is z' k = if the event at Zk is exogenous 



and 



-dx fk ^ 



ogk °Sk yl/ _ 
dd d. 



+ (29) 



if the event at Zk is endogenous (i.e., gk (x(6,Zk) ,6) = 0) and defined as long as ^fk(^k ) ^ 0- 
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IPA equations. To clarify the presentation, we first note that i = 1, ... ,M is used to index the points 
where uncertainty is measured; E, = 1, . . . ,T indexes the components of the parameter vector; and k = 
l,...,K indexes event times. In order to apply the three fundamental IPA equations ([27])-([29]> to our 
system, we use the state vector x(t) = [s(t) ,R\(t),. . . ,RM(t)] T and parameter vector 6 = [0\,..., dr] T - 
We then identify all events that can occur in Fig. [2] and consider intervals [t^(9),t^+i(9)) over which 
the system is in one of the 14 states shown for each i = 1, . . . ,M. Applying (27 1 to s(t) with f k (?) = 1 or 
-1 due to |lj and (Jl2j), the solution yields the gradient vector Vs(t) = [j^(t), ■ ■ ■ , jg^(0] T > where 



ds 



90. 



-(t^ + ), for? G [z k ,z k 



+ i. 



(30) 



for all k = 1, . . . ,K, i.e., for all states z(?) G {1, . . . , 14}. Similarly, let VRt(t) = [f§(0, • • • , |§^0)] T 
for i = 1,. . . ,M. We note from ^ that f k (?) = for states z(t) G Z\ = {7, 14}; f k (?) = A t for states 
z(t) G Z2 = {1,6,8, 13}; and f k (?) = A, — Bpi{s(t)) for all other states which we further classify into 
Z3 = {2,3, 11, 12} and Z4 = {4,5,9, 10}. Thus, solving (27i and using (30 1 gives: 



dRj 



VRi(t)=VRi(T2 



ifz(?)GZiUZ 2 
B ( d ^ W s (t+) • (t - Z k ) otherwise 



where dp '^ 

ds 

state. 



±i as evaluated from (4 ) depending on the sign of a,- — s(t) at each associated automaton 



We now turn our attention to the determination of Vs f T, H 



and VRi(xi~) which are needed to evaluate 
VRi (?) above. To do so, we use (28 1, which involves the event time gradient vectors Vz k = [j^,- ■ • , fgp] T 



for k = l,...,K (the value of K depends on T .) Looking at Fig. |2j there are three readily distinguishable 
cases regarding the events that cause discrete state transitions: 

Case 1: An event at time z k which is neither /?, = nor s = 8^, for any ^ = 1, . . . ,T. In this case, it is 
easy to see that the dynamics of both s(t) and /?,(?) are continuous, so that f k _i(%7) = f k (z^) in (28 1 
applied to s (?) and Ri(t), i = 1, . . . ,M gives: 



Vs(z; 



VR i (z+) = VR i (z k -), i=l, 



M 



(31) 



Case 2: An event /?; = at time z k . This corresponds to transitions 3 — > 7, 4 — > 7, 10 — > 14 and 11 — > 
14 in Fig. [2] where the dynamics of s(t) are still continuous, but the dynamics of Ri(t) switch from 
fk-\{\) = A, - Bpi{s(\)) to fk{*k) = ®- Thus, Vs (z k ) = Vs (z£) , but we need to evaluate z' k to 
determine V/?;(t^~). Observing that this event is endogenous, d29|) applies with g k = Rj = and we get 



dz k del ( X k ) 

£ = .,T,k = l,...,K 



90^ Ai-Bpiis^))' 



It follows from (28 (that 



,+sBRl, [A ; -flA(^-))]j|(T fe -) _ 

dd^ Tk) de^ Tk) Ai-Bp^)) 
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Thus, whenever an event occurs at %k such that Ri(ik) becomes zero, ^ (t^) is always reset to 
regardless of f| (t~). 

Case 3: An event at time Zk due to a control sign change at s = 0£, £ = 1, . . . ,r. This corresponds to 
any transition between the upper and lower part of the hybrid automaton in Fig. [2] In this case, the 
dynamics of Ri(t) are continuous and we have (t^ ) = (t^t) for all i,%,k. On the other hand, we 



have s(t k ) = u(z k ) = —u{x k ) = ±1. Observing that any such event is endogenous, (29l applies with 
gk = s— d?- = for some £ = 1 , . . . , T and we get 



I for) 



(32) 



Combining (32 1 with (28 1 and recalling that ^(t^) 



-ui/l, ), we have 



1 



d6? 



U T,. 



where ^ ) = because (0) = = (?) for all ? G [0, T,t), since the position of the agent cannot 
be affected by prior to this event. 



In this case, we also need to consider the effect of perturbations to 6j for j < i.e., prior to the current 
event time %k (clearly, for j > = since the current position of the agent cannot be affected 

by future events.) Observe that since gk = s — d^ = 0, we have ^ = for j / ^ and (29 1 gives 



ddj 



d, 
de 



so that using this in ( 28 1 we get 



ds + ds 



[«{*t) to _ ds 



de 



(O 



Combining the above results, the components of Vs(r^) where Zk is the event time when s(Zk) = 0£ for 
some ^ , are given by 



2 




if J = 1, 

if 7 = 5 
ifj>^ 



(33) 



It follows from (30 1 and the analysis of all three cases above that jg- (?) for all £, is constant throughout 



an optimal trajectory except at transitions caused by control switching locations (Case 3). In particular, 
for the Mi event corresponding to s(Tk) = 0?, t £ [t^, T], if u (?) = 1, then 4§- (?) = — 2 if ^ is odd, and 



(?) = 2 if E, is even; similarly, if m (?) 



•l,then^(?) 



2 if <^ is odd and jg- (?) = -2 if £ is even. 



In summary, we can write: 



dd> 



(0 



■lf-2u (?) 




t > ? k 
t < %k 



(34) 
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Finally, we can combine ( j34| with our results for |^ (?) in all three cases above. Letting s(t/) = 0^, we 



obtain the following expression for ||i (?) for all & > / , ? G [t^, Tk+i): 



507 



^7 W + < 



o ifz(0eziuz 2 

-l)« + 1 f M (T+).(?-T,) if z(?)ez 3 
^ -(-l)« +1 f M (T+)-(r-T,) ifz(?)eZ 4 



(35) 



with boundary condition 



ifz(T+)GZ! 

f£(T t -) otherwise 



(36) 



Objective Function Gradient Evaluation. Based on our analysis, the objective function ([7]) in problem 
PI can now be written as 7(0), a function of instead of u (?) and we can rewrite it as 



1 M K f%+l {6) 

^) = T II Ri(t,e) 



J? 



where we have explicitly indicated the dependence on 0. We then obtain: 



i M £ 



V/?; (?)rf? +Ri (Dfc+l) VT fc+ i (T ft ) Vt* 



Observing the cancelation of all terms of the form /?, (t*) Vt* for all k (with To = 0, Tk+i = T fixed), we 
finally get 

1 M K f r M (8) 

1 i= l k= QJ-Ct(ff) 



)dt. 



(37) 



The evaluation of V/(0) therefore depends entirely on V7? ( (?), which is obtained from (35 1-(|36|> and the 
event times z^, k = 1, . . .,K, given initial conditions s (0) = 0, Ri (0) for ? = 1, . . . ,M and V/?,-(0) = 0. 
Since V/?,- (?) itself depends only on the event times k = l,...,K, the gradient V/(0) is obtained by 
observing the switching times in a trajectory over [0, T] characterized by the vector 0. 



3.2.2 Multi agent solution where a > and b < L 

Next, we extend the results obtained in the previous section to the general multi-agent problem where 
we also allow a > and b < L. Recall that we require < a < r„ and L — r m <b< L, for at least 
some n,m= l,...,N since, otherwise, controlling agent movement cannot affect /?;(?) for all a,- located 
outside the sensing range of agents. We now include both parameter vectors 0„ = [0 n ,i, 0n,r„] T and 
w n = [w n .i, ■ ■ ■ Wn,r„V f° r eacn a g ent n and, for notational simplicity, concatenate them to construct = 
[01 , . . . , 0^v] T an d w = \w\ , . . . , %] T . The solution of problem PI reduces to the determination of optimal 
parameter vectors 6* and w* and we will use IPA to evaluate V/(0,w) = [ dJ ^^ dJ ^' w ^ ] J '. Similar to 

(pj, it is clear that this depends on Vfl/(f) = ' 
on a trajectory over [0, T] with given and w. 



99 dv 



and the event times Tj c ,k=l,...,K, observed 
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IPA equations. We begin by recalling the dynamics of /?, (t)in^ which depend on the relative positions 
of all agents with respect to a,- and change at time instants % k such that either Ri(xt) = with #,-(Tr) > 
or A, > BPi (s(Tk)) with/?,(T^) = 0. Moreover, using (TTb and our analysis in Section 3.1 the dynamics of 



s n (t), n=l,... ,N, in an optimal trajectory can be expressed as follows. Define ©„* = (0 n n f) if 
t, is odd and 0„ % = (6 n p , 6 n ,E,-\) if % i s even to De tne £th interval between successive switching points 
for any n = l,...,N, where d n Q = s n (0). Then, for E, = 1,2, . . ., 



Sn (t) 



1 s n (t) 6 8„ )? , £ odd 
-1 *„(/)€ 0^, £ even 
otherwise 



(38) 



where transitions for s n (?) from ±1 to =pl are incorporated by treating them as cases where w n p = 0, 
i.e., no dwelling at a switching point Q n p (in which case s n (?) = 0.) We can now concentrate on all 
events causing switches either in the dynamics of any /?,•(?),/ = 1, ... ,M, or the dynamics of any s n {t), 



n = 1, . . . ,/V. From (28 ), any other event at some time z k in this hybrid system cannot modify the values 

.... 



of VRi{t) 



\dSi(t) 




T 

or Vs„(t) = 


'ds„(t) ds n (t)~ 


de 


dw 


d 9 n dw n 



at t = T k . 



First, applying (27 ) to s n (t) with (t) = 1, — 1 or due to (|38j), the solution yields 

Vs„(t) = Vs„(tf), for ? G [lk^k+\) 



(39) 



for all k = 1 , . . . , K, n 

dRi , x dR 



I, 







de n £ 

and 

dRi 

dw n £ 



de, 



N. Similarly, applying (27 1 to Ri (?) and using ([6]) gives: 

ifRi(t)=0, Ai<BPi(s(t)) 
otherwise 



(0 



dRi 

dw n ^ 







(40) 

if/?,(?) = 0, Ai<BPi(s{t)) 
otherwise 

(41) 



Thus, it remains to determine the components of Vs„ (x^ ) and V/?/(t^~) in ( 39 1-(41 1 using (28 1. This 



involves the event time gradient vectors Vz k = for k = 1, . . . ,K, which will be determined 

through ( f29] ). There are three possible cases regarding the events that cause switches in the dynamics of 
Ri (?) or s„(t) as mentioned above: 

Case 1: An event at time %^ such that 7?,- (?) switches from /?; (?) = to Ri (?) = A,- — BPi (s(?)). In this 
case, it is easy to see that the dynamics of both s n (t) and Ri{t) are continuous, so that fk-\{\) = fk{^ k ) 



in (28 1 applied to s n (?) and Ri(t), i = 1, . . . ,M, n = 1, . . . ,/V, and we get 

= 1, 



Vj b (t+ 



V/? ! -(T, + )=V/? ! -(Tr) ! i=l ! 



.,/V 



(42) 
(43) 



Case 2: An event at time T k such that/?,- (f) switches from/?, (?) =Aj — BPj (s(?)) to/?, (?) = 0, i.e., /? ; -(t^) 
becomes zero. In this case, we need to first evaluate Vt^ from (29 1 in order to determine VR^tt) through 
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(28 1. Observing that this event is endogenous, (29 1 applies with gk = Ri = and we get 



(44) 



It follows from (28 (that 



V/?,-(T+)=V/? f (T i -)- 



[A,-(T t -)-^-(s(0)]v/e,-(T t -) 

A ; .(t-)-^(t-) 



(45) 



Thus, VRj(r£) is always reset to regardless of VRj(r k ). In addition, (|42j) holds, since the the dynamics 
of s n (t) are continuous at time % 

Case 3: An event at time such that the dynamics of s n {t) switch from ±1 to 0, or from to dbl. 
Clearly, ( pOj ) holds since the the dynamics of R,-(t) are continuous at this time. However, determining 

ds fr + ) 

Vs„ (t^J is more elaborate and requires us to consider its components separately, first "^ g k ' and then 

dw„ 



Case 3.1: Evaluation of 



MO 

de„ ■ 



Case 3.1.1: An event at time such that the dynamics of s n (t) in (38 1 switch from ±1 to 0. This is an 
endogenous event and ( 29 ) applies with gk = s n — 6 n ^ = for some £ = 1 , . . . , T n and we have: 



de, 



and (P28l yields 



ds n , — I — \ 



ds n 



1 



) + [««(% )-°]- 



" 3^ ) 



(46) 



(47) 



As in Case 3 of Section 3.2.1 we also need to consider the effect of perturbations to 6j for j < i.e., 
prior to the current event time (clearly, for 7 > = since the current position of the agent 

cannot be affected by future events.) Observe that = 0, therefore, (29 ) becomes 



and using this in (28 ) gives: 

ds n 



de nJ 



dd nJ 



n i - 



de, 



(48) 







(49) 



Thus, combining the above results, when s q (Zk) 
we have 

ds n . — I — \ 



B„ i for some E, and the agent switches from ±1 to 0, 

0, ifjv« (50) 

1, if / = § (5U) 
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Case 3.1.2: An event at time Zk such that the dynamics of s n (t) in (38 1 switch from to ±1. This is 



an induced event since it is triggered by the occurrence of some other endogenous event when the agent 
switches from ±1 to (see Case 3.1.1 above.) Suppose the agent starts from an initial position s n (0) = a 
with u n (0) = 1 and T& is the time the agent switches from the to ±1 at the switching point 6 n t . If B n t 



is such that u n \z£ 



1, then ^ is even and z^ can be calculated as follows: 



Tjfc = -a)+W n s + (6 n ,\ -Onl) +W n ,2 + --- + {^n,t,-\ ~ + w n,£, 

( £-1 S-2 \ % 



,v=l, v odd 



v=2, v even 



v=l 



(51) 



Similarly, if n £ is the switching point such that u„ (t^) = — 1, then E, is odd and we get: 



^ = 2 £ e„ )V - £ e„ ;V J + J^w^v+e^ 

\v=l,vodd v=2.veven / v=l 



We can then directly obtain -J^- as 



(52) 



50, 



-J«w(m(t^)) 



(53) 



Using (53 1 in (28 1 gives: 



ds n 



dd, 



(z k ) + [0-u{z+)]-[-sgn(u(z+))} 



ds„ 



(54) 



Once again, we need to consider the effect of perturbations to 6j for j < i.e., prior to the current event 
time Zk (clearly, for j > = 0.) In this case, from (51 >-(52>, we have 



^t=2, if ; odd 
4^- = -2, if /even 



(55) 



and it follows from ( 28 1 that for j < t, : 



dOnJ 



_d_ 

de, 
d_ 
de. 



odd 



^ (z k ) + 2, if k„ (t^) = 1 , 7 even, or m„ (t+) = -1 , j 
, g-(Tj~)-2, if «„ (t^ + ) = 1 , j odd, or u n (t+) = - 1 , ; even 



(56) 



Case 3.2: Evaluation of 



Case 3.2.1: An event at time such that the dynamics of s n (t) in (38 1 switch from ±1 to 0. This is an 



endogenous event and (29 1 applies with gk = s n — 6 n ^ = for some t, = 1, . . . ,r„. Then, for any j < 
we have: 



dz k 



W n ,j ) 
d\V n j U n (z^) 



(57) 
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Combining (57 1 with (28 1 and since u n (% k ) = ±1, we have 
ds„ 



1W, 



n.j 



dw n ,, 



(z k ) + [u n (z k 



0] 



(58) 



J.i.2, T/t is given by ( 



Using this result in ( 28 



Case 3.2.2: An event at time such that the dynamics of s n {t) in (38 ) switch from to ±1. As in Case 
( 5 1 1 or ( 52 1, depending on the sign of u q (t^) . Thus, we have ^— - = 1, for j < % . 



ds n 

dw n .j 



and observing that 5^*7 Ofy ) = from (58 1, we have 
ds 



i XV 



^(t-) + [0- I4 „(t+)].1 



"■j 



< »forj<^ 



(59) 



Combining the above results, we have for Case 3.2: 



IW n . 



0, if u„ (x k ) = ±1, u„ (t^) =0 
=Fl, if u„ (i k ) = 0, u„ (t+) = ±1 



(60) 



Finally, note that g^r-(0 = for t £ [0, Tjt), since the position of the agent n cannot be affected by w n £ 
prior to such an event. 



Objective Function Gradient Evaluation. Proceeding as in the evaluation of V/(0) in Section 3.2.1 
we are now interested in minimizing the objective function J(d, w) in ([7]) with respect to 6 and w and we 
can obtain VJ(d,w) 



\dJ(9,w) dj{6.w)^j 



de 



dw 



as 



M K 



VJ(0,w) 



T*+l(0,w) 



££/ v Ri (t)dt 



This depends entirely on VR,-(t), which is obtained from d40b and (41 1 and the event times T^, k 



1 , . . . , K, given initial conditions s n (0) = a for n = 1 , . . . , Af, and (0) for i = 1 , . . . , M. In ( 40 



ds„ ( T + 

is obtained through (|43|> and (|45|>, whereas k ' is obtained through (|39|>, (1421), (|50|>, and ( 



(t^) is again obtained through (43 1 and (45 1, whereas 3s ^ k - is obtained through (|42|), and ((60 



Ml ( z +) 



561. In (41 1, 



Remark 2. Observe that the evaluation of V7?, (?), hence VJ(6,w), is independent of A,, i = 1, . . . ,M, 
i.e., the values in our uncertainty model. In fact, the dependence of VRi (t) on A,-, / = 1 , . . . ,M, manifests 
itself through the event times Xk, k = 1, . . . ,K, that do affect this evaluation, but they, unlike A; which 
may be unknown, are directly observable during the gradient evaluation process. Thus, the IPA approach 
possesses an inherent robustness property: there is no need to explicitly model how uncertainty affects 
Ri(t) in Consequently, we may treat A, as unknown without affecting the solution approach (the 
values of VRi (t) are obviously affected). We may also allow this uncertainty to be modeled through 



random processes {A(?)}, J = 1, - ■ ■ ,M; in this case, however, the result of Proposition 3.3 no longer 



applies without some conditions on the statistical characteristics of {A,-(?)} and the resulting V/(0, w) is 
an estimate of a stochastic gradient.) 
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3.3 Objective Function Optimization 



We now seek to obtain 6* and w* minimizing J(6,w) through a standard gradient-based optimization 
scheme of the form 

[e'+V +1 ] T = [eV] 1 -^ r) w ]vj{e l ,w l ) (6i) 

where {Hwl ^ appropriate step size sequences and VJ(d' ,w') is the projection of the gradient 

VJ(8 1 ,w ! ) onto the feasible set (the set of 6 I+l satisfying the constraint p6] >, a < 8 1+l < b, and w > 0). 
The optimization scheme terminates when |V/(0,w)| < £ (for a fixed threshold e) for some 6 and w. 
Our IPA-based algorithm to obtain 6* and w* minimizing J(6, w) is summarized in Algorithm [I] where 
we have adopted the Armijo method in step-size selection (see |23]) for { [rj' g tj/J }. 



One of the unusual features in (61 1 is the fact that the dimension T* of 6* and w* is a priori unknown (it 
depends on T). Thus, the algorithm must implicitly determine this value along with 6* and w*. One can 
search over feasible values of Y n € {1,2, . . .} by starting either with a lower bound Y n = 1 or an upper 
bound to be found. The latter approach results in much faster execution and is followed in Algorithm 
[T] An upper bound is determined by observing that d n £ is the switching point where agent n changes 
speed from 1 to for E, odd and from — 1 to for E, even. By setting these two groups of switching 
points so that their distance is sufficiently small and waiting times w n = for each agent, we determine 
an approximate upper bound for r„ as follows. First, we divide the feasible space [a, b] evenly into N 
intervals: [a + ^ (b — a) ,a + S (b — a)], n = l,...,N. Define D n = a+ (b — a) to be the geometric 
center of each interval and set 

=D n -o if t, even 
Q nA =D n + o if § odd K } 

so that the distance between switching points 6 n t for odd and even is 2a, where a > is an arbitrarily 
small number, n = 1, . .. ,N. In addition, set w n = 0. Then, T must satisfy 



0»,i -s„ (0) + 2a (r„ -1) < T < Qn,x-sn (0) +2oT n (63) 
n = l,...,N, where F n is the number of switching points agent n can reach during (0, T], given 6 n £ are 



defined in ( |62[ >. From ( |63| > and noting that r„ is an integer, we have 

1 



2a 



[T-^+SntO)] 



(64) 



where [•] is the ceiling function. Clearly, reducing a increases the initial number of switching points T n 
assigned to agent n and T n — > oo as a — > 0. Therefore, a is selected sufficiently small while ensuring that 
the algorithm can be executed sufficiently fast. 

As Algorithm [I] repeats steps 3-6, w n £ > and distances between B n ^ for E, odd and even generally 
increase, so that the number of switching points agent n can actually reach within T decreases. In other 
words, as long as a is sufficiently small (hence, T n is sufficiently large), when the algorithm converges 
to a local minimum and stops, there exists < Y n , such that 6 n r is the last switching point agent n 
can reach within (0, T], n = 1, . .. ,N. Observe that there generally exist such that < E, < F n which 
correspond to points 6 n % that agent n cannot reach within (0, T]; the associated derivatives of the cost 
with respect to such 8 n % are 0, since perturbations to these Q n * will not affect s„ (t), t G (0, T] and 
thus the cost J(6,w). When |V/(0,w)| < e, we achieve a local minimum and stop, at which point the 
dimension of 6* and w* is £„. 



21 



Algorithm 1 : IPA-based optimization algorithm to find 6* and w* 



l: Pick a > and e > 0. 

2: Define D n =a + (b — a),n = l,...,N, and set 



* = D n — a if § even 
re | =D n + o if ^ odd 



Set w = [w b . . . ,w N ] = 0., where w„ = [w n> i,. . . ,w„^J and T„ = [7- 6 nA +s n (0)]] 
repeat 

Compute s n (t), t G [0,T] using s„(0), (12 1, and w for n = 1, 
Compute V/(0,w) and update 0,w through ( |6Tj ) 
until |V7(0,w)| < £ 



7: Set0* 



and w* 



switching point agent n can reach within (0, T],n- 



, where is the index of 6„ t, which is the last 
l,...,N 



4 Numerical Results 



In this section we present some examples of persistent monitoring problems in which agent trajectories 
are determined using Algorithm[T] The first four are one-agent examples with L = 20, M = 21, (X\ = 0, 
(Xm = 20, and the remaining sampling points are evenly spaced over [0,20]. The sensing range in ^ 
is set to r = 4, the initial values of the uncertainty functions in ^ are /?,(0) = 4 for all i, and the time 
horizon is T = 400. In Fig. [3ja) we show results where the agent is allowed to move over the entire space 
[0,20] and the uncertainty model is selected so that B = 3 and A, = 0.1 for all i = 1, . . . ,20, whereas in 
Fig. [5Jb) the feasible space is limited to [a, b] with a = r = 4 and b = L — r = 16. The top plot in 
each example shows the optimal trajectory s*{t) obtained, while the bottom shows the cost J(d l ,w 1 ) as 
a function of iteration number. In Fig. [4j the trajectories in Fig. [3ja),(b) are magnified for the interval 
t G [0,75] to emphasize the presence of strictly positive waiting times at the switching points. In addition, 
maximum, minimum and mean values for the uncertainty function of each sampling point in these two 
cases are shown in Fig. [5] Observe that when a = r and b = L — r, an instability arises at the last two 
sampling points of both ends of the mission space; this is expected since the agent's sensing range can 
only marginally reach the two end points from a = 4 and b = 16. 

In Fig. [3jc) we show results for a case similar to Fig. [3ja) except that the values of A; are selected so 
that Ao = A20 = 0.5, while A,- = 0. 1 , i = 1 , . . . , 19. We should point out that even though it seems that the 
trajectory includes switching points at the two end points, this is not the case: the switching points are 
very close but not equal to these end points, consistent with Proposition 3.1 In Fig. |3jd), on the other 



hand, the values of A, are allowed to be random, thus dealing with a persistent monitoring problem in a 
stochastic mission space. In particular, each A is treated as a piecewise constant random process (A,(f)} 
such that Aj(t) takes on a fixed value uniformly distributed over (0.075,0.125) for an exponentially 
distributed time interval with mean 10 before switching to a new value. Note that the behavior of the 
system in this case is very similar to Fig. [3ja) where A, = 0.1 for all i = 1, ... ,20 without any change 



in the way in which V/(0',w') is evaluated in executing (61 1. As already pointed out, this exploits a 



robustness property of IPA which makes it independent of the values of A, . In general, however, when 



Ai(t) is not time-invariant, Proposition 3.3 may no longer apply, since an extra term A (?) would be 
present in (22i. In such a case, u* (t) may be nonzero when A* (t) = and the determination of an 
optimal trajectory through switching points and waiting times alone may no longer be possible. In the 
case of [5Jd), A,(?) changes sufficiently slowly to maintain the validity of Proposition 3.3 over relatively 
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long time intervals, under the assumption that w.p. 1 no event time coincides with the jump times in any 

In all cases, we initialize the algorithm with a = 5 and e = 2 x 10~ 10 . The algorithm running times are 
approximately 10 sec using Armijo step-sizes. Note that although the number of iterations for the exam- 
ples shown may substantially vary, the actual algorithm running times do not. This is simply because the 
Armijo step-size method may involve several trials per iteration to adjust the step-size in order to achieve 
an adequate decrease in cost. In Fig. |3ja),(d), red line shows J vs. number of iterations using constant 
step size and they almost converges to the same optimal value. Non-smoothness in Fig. [3jd) comes from 
the fact that it is a stochastic system. Note that in all cases the initial cost is significantly reduced indicat- 
ing the importance of optimally selecting the values of the switching points and associated waiting times 
(if any). 

Figure [6] shows two two-agent examples with L = 40, M = 41 and evenly spaced sampling points over 
[0,L], A, ■ = 0.01 , B = 3, r = 4, # ; (0) = 4 for all i and T = 400. In Fig. ga) the agents are allowed to move 
over the whole mission space [0,L], while in Fig. |6jb) they are only allowed to move over [a, b] where 
a = r and b = L — r. We initialize the algorithm with the same a and e as before. The algorithm running 
time is approximately 15 sec using Armijo step-sizes, and we observe once again significant reductions 
in cost. 



5 Conclusion 
References 

[1] I. Rekleitis, V. Lee-Shue, A. New, and H. Choset, "Limited communication, multi-robot team based 
coverage," in Robotics and Automation, 2004. Proceedings. ICRA'04. 2004 IEEE International 
Conference on, vol. 4. IEEE, 2004, pp. 3462-3468. 

[2] J. Cortes, S. Martinez, T. Karatas, and F. Bullo, "Coverage control for mobile sensing networks," 
Robotics and Automation, IEEE Transactions on, vol. 20, no. 2, pp. 243-255, 2004. 

[3] W. Li and C. Cassandras, "A cooperative receding horizon controller for multivehicle uncertain 
environments," IEEE Transactions on Automatic Control, vol. 51, no. 2, pp. 242-257, 2006. 

[4] A. Girard, A. Howell, and J. Hedrick, "Border patrol and surveillance missions using multiple 
unmanned air vehicles," in 43rd IEEE Conference on Decision and Control, vol. 1. IEEE, 2005, 
pp. 620-625. 

[5] B. Grocholsky, J. Keller, V. Kumar, and G. Pappas, "Cooperative air and ground surveillance," 
IEEE Robotics & Automation Magazine, vol. 13, no. 3, pp. 16-25, 2006. 

[6] S. Smith, M. Schwager, and D. Rus, "Persistent robotic tasks: Monitoring and sweeping in changing 
enviroments," IEEE Transactions on Robotics, 2012, to appear. 

[7] D. Paley, F. Zhang, and N. Leonard, "Cooperative control for ocean sampling: The glider coor- 
dinated control system," IEEE Transactions on Control Systems Technology, vol. 16, no. 4, pp. 
735-744,2008. 



23 



20 



Agent position vs. time 



Agent position vs. time 




50 100 150 200 250 300 350 400 
Timet 

Cost J vs. Number of iterations 



100 



8 50 
o 




- Armijo step 

- costant step 



10 20 30 40 50 60 
Number of iterations 



70 



(a) a = 0, b = 20. A t = 0.1, i = 1, . . . , 20. /* = 17.77. 

Agent position vs. time 




50 1 00 1 50 200 250 300 350 400 
Time t 

Cost J vs. Number of iterations 



6 8 10 12 14 16 18 
Number of iterations 




50 1 00 1 50 200 250 300 350 400 
Time t 

Cost J vs. Number of iterations 




50 1 00 1 50 200 250 300 350 400 
Time t 

Cost J vs. Number of iterations 



60 
50 
2 40 
8 30 
20 
10 




20 30 40 50 

Number of iterations 



(c) a = 0,* = 20. A = A 20 = 0.5,A,- = 0.1, i = 1, . . . , 19. (d) a = 0,b = 20. A,- (Af,) - V (0.075,0.125), At t - O.le 
J* = 39.30. J* = 17.54. 



-0.1/ 



Figure 3: One agent example. L = 20, T = 400. For each example, top plot: optimal trajectory; bottom 
plot: J versus iterations. 
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(a)a = 0,fo = 20. (b)a = 4,fc=16. 

Figure 4: Magnified trajectory for sub-figure (a) and (b) in Fig. [5J t G [0,75]. 
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(a)a = 0,fc = 20. (b) a = 4,6=16. 



Figure 5 : Max, min and mean uncertainty value for each sampling point. 
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(a) a = 0,b = 20. J* = 17.77. (b) a = 4,b = 16. J* = 39.14. 

Figure 6: Two agent example. L = 40, T = 400. Top plot: optimal trajectory. Bottom plot: J versus 
iterations. 
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